Since release of the Biological Opinion on oil and gas activities in the Gulf of Mexico (NMFS 2020) that used a published density surface model (Roberts et al. 2016) to describe the distribution of the critically endangered Rice’s whale (Balaenoptera ricei), a new density surface model (Litz et al. 2022) has been made available. Importantly, this model extends the distribution of Rice’s whale beyond its initial core habitat in the Eastern Gulf of Mexico to the West, where it had previously only been acoustically detected (Soldevilla et al. 2022). This report replicates the Biological Opinion’s ship strike analysis using the newer Rice’s whale distributional model. Given the wider distribution of Rice’s whale, an alternative new Whale Area is suggested to reduce ship strike risk with the Rice’s whale based simply on location (25.5º N and higher) and depth (100 to 400 m).
2 Whale Densities
The new density surface model (Litz et al. 2022) uses approximately 40 km2 hexagons as its spatial unit to describe number of individuals per 40 km2 in a Lambert Conformal Conic projection, whereas the original model (Roberts et al. 2016) used 100 km2 cells in a custom equal area Albers projection to describe number of individuals per 100 km2. The spatial unit for this new analysis is also 100 km2 cells but in the web Mercator projection (EPSG:3857) in order to readily map results online with common “slippy” basemaps, like the Esri Ocean Basemap. All layers were clipped to the study area of the U.S. Exclusive Economic Zone (EEZ) within the Gulf of Mexico.
Normally converting polygons to raster extracts only the centroid point of the raster cell from the underlying polygon. In order to capture the entirety of the underlying geometric densities, a vector-based intersection was first performed on all layers (whale hexagons, ship cells, and new units) before summarizing to the raster cell as area-weighted means.
In order to adjust for slight differences from projecting coordinate reference systems and rounding errors, the new 100 km2 whale density grid was adjusted so the sum of individuals predicted throughout the study area is equal to 51.3, the most recent abundance estimate (Garrison, Ortega-Ortiz, and Rappucci 2020) from the same 2017 and 2018 surveys as the new density model (Litz et al. 2022) was derived.
Compared to the original distribution (Figure 1), the whales are now concentrated along the strip from 100 to 400 m extending into the Western Gulf of Mexico (Figure 2).
Figure 1: Map of previous whale densities (Roberts et al. 2016) as 100 km2 cells used by (NMFS 2020) showing the dominance in the northeastern corner of the Gulf of Mexico. The original Whale Area (p. 292 of NMFS 2020) is depicted by the pink outline polygon for vessel slowdown and nighttime avoidance. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
The original Whale Area is described (p. 292 of NMFS 2020) as:
This opinion defines the Bryde’s whale area to include the area from 100- to 400- meter isobaths from 87.5° W to 27.5° N as described in the status review (Rosel 2016) plus an additional 10 km around that area.
There was only a tiny marginal improvement in capturing additional whale densities by adding the 10 km buffer used to create the original Whale Area (pink outline in Figure 1). In generating the new Whale Area, the ease of navigation with simpler description in terms only of a southern limit and depth range outweighed this marginal improvement (red outline in Figure 2). This new Whale Area captures 94% of the population from the new density estimates (Litz et al. 2022) compared to only 52% of the original Whale Area (Table 1).
Figure 2: Map of new whale densities (Litz et al. 2022) as 100 km2 cells showing a distribution throughout the region. The newly recommended Whale Area is depicted by the red outline polygon for vessel slowdown and nighttime avoidance using similar logic as to (NMFS 2020). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Table 1:
Table of new whale densities (Litz et al. 2022) summarized by total
study area (U.S. Gulf of Mexico), previous Whale Area (NMFS 2020) and
newly proposed Whale Area.
Item
#
%
Whales in Study (U.S. Gulf of Mexico)
51
100%
Whales in Original Whale Area (NMFS, 2020)
27
52%
Whales in New Whale Area
48
94%
3 Vessel Traffic
In order to evaluate the threat of ship strike to Rice’s whales, we used the same AIS data from 2014 to 2018 as the Biological Opinion (NMFS 2020). This data is based on a grid of cells ~126 km2 in Albers equal area projection. Traffic in terms of kilometers (km) traversed within a cell was differentiated based on speed (≤ 10 knots or > 10 knots) and type (oil & gas or all types).
Figure 3: Map of annual average traffic (km) for all vessel types at all speeds from AIS data (2014 to 2018).
Figure 4: Map of annual average traffic (km) for oil and gas vessels at all speeds from AIS data (2014 to 2018).
Figure 5: Map of annual average traffic (km) for all vessel types > 10 knots from AIS data (2014 to 2018).
Figure 6: Map of annual average traffic (km) for oil and gas vessels > 10 knots from AIS data (2014 to 2018).
4 Vessel Risk to Whales
Figure 7: Map of risk (# whales * km vessel traffic) for all vessels at all speeds.
Figure 8: Map of risk (# whales * km vessel traffic) for oil and gas vessels > 10 knots.
The overall vessel strike risk to Rice’s whale is presented in Table 2 [similar to Table 49 (NMFS 2020)].
Table 2:
Vessel strike risk (# whales * km vessel traffic) to Rice’s whales
for oil and gas vessels compared with all vessels.
Reduction of vessel strike risk from oil and gas vessel traffic (#
whales * km vessel traffic) to Rice’s whales with enforcement of
proposed mitigation areas for oil and gas vessels.
Year
Risk Reduction by Area
Original
%
New
%
All speeds - All vessels
2015
11,025
12%
86,016
93%
2016
10,963
13%
79,186
93%
2017
10,343
12%
80,252
93%
2018
11,956
12%
93,215
93%
Avg
11,072
12%
84,667
93%
All speeds - Oil & Gas vessels
2015
2,516
3%
34,977
38%
2016
1,827
2%
30,287
35%
2017
2,175
3%
31,073
36%
2018
1,520
2%
35,309
35%
Avg
2,010
2%
32,911
36%
> 10 knots - All vessels
2015
7,842
11%
66,483
93%
2016
8,194
12%
63,161
93%
2017
7,644
11%
63,856
93%
2018
9,105
11%
74,307
93%
Avg
8,196
11%
66,952
93%
> 10 knots - Oil & Gas vessels
2015
949
1%
21,089
29%
2016
659
1%
19,655
29%
2017
721
1%
20,145
29%
2018
450
1%
23,116
29%
Avg
695
1%
21,001
29%
5 Reproducible Results
References
Garrison, Lance, Joel Ortega-Ortiz, and Gina Rappucci. 2020. “Abundance of Marine Mammals in Waters of the U.S. Gulf of Mexico During the Summers of 2017 and 2018.” PRBD-2020-07.
Litz, Jenny, Laura Aichinger Dias, Gina Rappucci, Anthony Martinez, Melissa Soldevilla, Lance Garrison, and Keith Mullin. 2022. “Cetacean and Sea Turtle Spatial Density Model Outputs from Visual Observations Using Line-Transect Survey Methods Aboard NOAA Vessel and Aircraft Platforms in the Gulf of Mexico.”
NMFS. 2020. “‘Biological Opinion on the Federally Regulated Oil and Gas Program Activities in the Gulf of Mexico,’ 13 March 2020, a Consultation Conducted by the Endangered Species Act Interagency Cooperation Division, Office of Protected Resources, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, U.S. Department of Commerce.”
Roberts, Jason J., Benjamin D. Best, Laura Mannocci, Ei Fujioka, Patrick N. Halpin, Debra L. Palka, Lance P. Garrison, et al. 2016. “Habitat-Based Cetacean Density Models for the U.S. Atlantic and Gulf of Mexico.”Scientific Reports 6 (March): 22615. https://doi.org/10.1038/srep22615.
Soldevilla, Melissa S., Amanda J. Debich, Lance P. Garrison, John A. Hildebrand, and Sean M. Wiggins. 2022. “Rice’s Whales in the Northwestern Gulf of Mexico: Call Variation and Occurrence Beyond the Known Core Habitat.”Endangered Species Research 48 (July): 155–74. https://doi.org/10.3354/esr01196.
Source Code
---title: "Spatial analysis of ship strike risk for Rice’s whale in the Gulf of Mexico"author: "Benjamin D. Best, Ph.D. (<ben@ecoquants.com>)"date: nowdate-format: "YYYY-MM-DD HH:mm (z)"bibliography: "ricei.bib"format: html: toc: true number-sections: true number-depth: 3 code-fold: true code-tools: true docx: toc: true toc-depth: 3 toc-title: "Contents" number-sections: true code-annotations: false execute: echo: false warning: falseeditor_options: chunk_output_type: console---```{r setup, include=FALSE}knitr::opts_chunk$set(echo = F,warning = F,message = F)source(here::here("scripts/functions.R"))source(here("scripts/paths.R"))ply_study <-read_sf(study_gpkg)ply_cells <-read_sf(cells_geo)ply_units <-read_sf(units_gpkg)ply_ships <-read_sf(ships_geo)ply_whales <-read_sf(whales_geo)ply_bia <-read_sf(bia_geo)ply_wab <-read_sf(wab_geo)ply_wan <-read_sf(wan_geo)tbl_ships <-read_csv(ships_csv)tbl_whales <-read_csv(whales_csv)lns_depth_contours <-read_sf(depth_contours_geo)stk_rast <-rast(rast_tif)r_cells <- stk_rast["cell_id"]tbl_units <- ply_units |>st_drop_geometry()```## AbstractSince release of the Biological Opinion on oil and gas activities in the Gulf of Mexico [@nmfsBiologicalOpinionFederally2020] that used a published density surface model [@robertsHabitatbasedCetaceanDensity2016] to describe the distribution of the critically endangered Rice's whale (_Balaenoptera ricei_), a new density surface model [@litzCetaceanSeaTurtle2022] has been made available. Importantly, this model extends the distribution of Rice's whale beyond its initial core habitat in the Eastern Gulf of Mexico to the West, where it had previously only been acoustically detected [@soldevillaRiceWhalesNorthwestern2022]. This report replicates the Biological Opinion's ship strike analysis using the newer Rice's whale distributional model. Given the wider distribution of Rice's whale, an alternative new Whale Area is suggested to reduce ship strike risk with the Rice's whale based simply on location (25.5º N and higher) and depth (100 to 400 m).## Whale DensitiesThe new density surface model [@litzCetaceanSeaTurtle2022] uses approximately 40 km^2^ hexagons as its spatial unit to describe number of individuals per 40 km^2^ in a Lambert Conformal Conic projection, whereas the original model [@robertsHabitatbasedCetaceanDensity2016] used 100 km^2^ cells in a custom equal area Albers projection to describe number of individuals per 100 km^2^. The spatial unit for this new analysis is also 100 km^2^ cells but in the web Mercator projection (EPSG:3857) in order to readily map results online with common "slippy" basemaps, like the Esri Ocean Basemap. All layers were clipped to the study area of the U.S. Exclusive Economic Zone (EEZ) within the Gulf of Mexico.Normally converting polygons to raster extracts only the centroid point of the raster cell from the underlying polygon. In order to capture the entirety of the underlying geometric densities, a vector-based intersection was first performed on all layers (whale hexagons, ship cells, and new units) before summarizing to the raster cell as area-weighted means.In order to adjust for slight differences from projecting coordinate reference systems and rounding errors, the new 100 km^2^ whale density grid was adjusted so the sum of individuals predicted throughout the study area is equal to `r n_whales_Garrison2020`, the most recent abundance estimate [@garrisonAbundanceMarineMammals2020] from the same 2017 and 2018 surveys as the new density model [@litzCetaceanSeaTurtle2022] was derived. <!-- - discuss changing distribution -->Compared to the original distribution (@fig-map-whales-old), the whales are now concentrated along the strip from 100 to 400 m extending into the Western Gulf of Mexico (@fig-map-whales-new).```{r}#| label: fig-map-whales-old#| fig-cap: "Map of previous whale densities [@robertsHabitatbasedCetaceanDensity2016] as 100 km^2^ cells used by [@nmfsBiologicalOpinionFederally2020] showing the dominance in the northeastern corner of the Gulf of Mexico. The original Whale Area [p. 292 of @nmfsBiologicalOpinionFederally2020] is depicted by the pink outline polygon for vessel slowdown and nighttime avoidance. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast(r =rast(whales_roberts2016_img), legend_title ="Whales<br><small>(Roberts, 2016)<br>(# / 100 km<sup>2</sup>)</small>",group ="whales_per_100km2_Roberts2016",add_depth_contours =T,add_ply_wab = T)```The original Whale Area is described [p. 292 of @nmfsBiologicalOpinionFederally2020] as:> This opinion defines the Bryde’s whale area to include the area from 100- to 400- meter isobaths from 87.5° W to 27.5° N as described in the status review (Rosel 2016) plus an additional 10 km around that area.There was only a tiny marginal improvement in capturing additional whale densities by adding the 10 km buffer used to create the original Whale Area (pink outline in @fig-map-whales-old). In generating the new Whale Area, the ease of navigation with simpler description in terms only of a southern limit and depth range outweighed this marginal improvement (red outline in @fig-map-whales-new). This new Whale Area captures 94% of the population from the new density estimates [@litzCetaceanSeaTurtle2022] compared to only 52% of the original Whale Area (@tbl-whales-by-area).```{r}#| label: fig-map-whales-new#| fig-cap: "Map of new whale densities [@litzCetaceanSeaTurtle2022] as 100 km^2^ cells showing a distribution throughout the region. The newly recommended Whale Area is depicted by the red outline polygon for vessel slowdown and nighttime avoidance using similar logic as to [@nmfsBiologicalOpinionFederally2020]. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast(r = stk_rast["whales_avg_n_per_100km2"], legend_title ="Whales<br><small>(Litz, 2022)<br>(# / 100 km<sup>2</sup>)</small>",group ="whales_per_100km2_Litz2022",add_depth_contours=T,add_ply_wan = T)``````{r}#| label: tbl-whales-by-area#| tbl-cap: Table of new whale densities (Litz et al. 2022) summarized by total study area (U.S. Gulf of Mexico), previous Whale Area (NMFS 2020) and newly proposed Whale Area.# Manual docx: mv tbl outside figure table. Reference > Insert Caption. Add Table 1 and Insert Hyperlink to This Document "tbl-whales-by-area" to replace "?@tbl-whales-by-area".a <-read_csv(whales_n_by_area_csv)d <-tribble(~Item, ~`#`, ~`%`,"Whales in Study (U.S. Gulf of Mexico)", a$n_whales_gom, 1,"Whales in Original Whale Area (NMFS, 2020)", a$n_whales_wab, a$pct_whales_wab,"Whales in New Whale Area", a$n_whales_wan, a$pct_whales_wan)d |>gt() |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent(`%`, decimals =0) |>opt_stylize(style =6)```## Vessel TrafficIn order to evaluate the threat of ship strike to Rice's whales, we used the same AIS data from 2014 to 2018 as the Biological Opinion [@nmfsBiologicalOpinionFederally2020]. This data is based on a grid of cells ~126 km^2^ in Albers equal area projection. Traffic in terms of kilometers (km) traversed within a cell was differentiated based on speed (≤ 10 knots or > 10 knots) and type (oil & gas or all types).```{r}#| label: fig-ships-avg-all-gt01#| fig-cap: "Map of annual average traffic (km) for all vessel types at all speeds from AIS data (2014 to 2018)."map_rast_jenks(stk_rast["ships_avg_all_gt01"], "Traffic,<br><small>all types,<br>all speeds,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-boem-gt01#| fig-cap: "Map of annual average traffic (km) for oil and gas vessels at all speeds from AIS data (2014 to 2018)."map_rast_jenks(stk_rast["ships_avg_boem_gt01"], "Traffic,<br><small>oil & gas,<br>all speeds,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-all-gt10#| fig-cap: "Map of annual average traffic (km) for all vessel types > 10 knots from AIS data (2014 to 2018)."map_rast_jenks(stk_rast["ships_avg_all_gt01"], "Traffic,<br><small>all types,<br>> 10 knots,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-boem-gt10#| fig-cap: "Map of annual average traffic (km) for oil and gas vessels > 10 knots from AIS data (2014 to 2018)."map_rast_jenks(stk_rast["ships_avg_boem_gt01"], "Traffic,<br><small>oil & gas,<br>> 10 knots,<br>avg year<br>(km)</small>", "RdYlGn")```## Vessel Risk to Whales```{r}#| label: fig-risk-avg-all-gt01#| fig-cap: "Map of risk (# whales * km vessel traffic) for all vessels at all speeds."stk_rast["risk_avg_all_gt01"] |>map_rast("Risk,<br><small>all types,<br>all speeds<br>(# * km)</small>", "RdYlGn")``````{r}#| label: fig-risk-avg-boem-gt10#| fig-cap: "Map of risk (# whales * km vessel traffic) for oil and gas vessels > 10 knots."stk_rast["risk_avg_boem_gt10"] |>map_rast("Risk,<br><small>oil & gas,<br>> 10 knots <br>(# * km)</small>", "RdYlGn")```The overall vessel strike risk to Rice's whale is presented in @tbl-risk-overview \[similar to Table 49 [@nmfsBiologicalOpinionFederally2020]\].```{r}#| label: tbl-risk-overview#| tbl-cap: "Vessel strike risk (# whales * km vessel traffic) to Rice's whales for oil and gas vessels compared with all vessels."d <-read_csv(risk_overview_csv) |>pivot_longer(-yr,names_to ="var",values_to ="val") |># table(d$var)separate_wider_delim( var,delim ="_",names =c("v","type","speed")) |>pivot_wider(names_from =c(v, type),names_sep ="_",values_from = val) |>relocate(speed) |>arrange(speed, yr)d <- d |>mutate(yr =as.character(yr)) |>bind_rows(bind_cols(yr ="Avg", d |>group_by(speed) |>summarise(across(where(is.numeric) &!yr, mean)) ) ) |>arrange(speed, yr)d <- d |>mutate(speed =case_match( speed, "gt01"~"All speeds","gt10"~"> 10 knots") )d |>group_by(speed) |>gt() |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent( pct_boem, decimals =0) |>opt_stylize(style =6) |>cols_label(yr ="Year",risk_all ="All Vessels",risk_boem ="Oil & Gas",pct_boem ="%") |>tab_spanner(label ="Vessel Strike Risk",columns =c( risk_all, risk_boem, pct_boem))```Reference @tbl-risk-reduction-by-areas.```{r}#| label: tbl-risk-reduction-by-areas#| tbl-cap: "Reduction of vessel strike risk from oil and gas vessel traffic (# whales * km vessel traffic) to Rice's whales with enforcement of proposed mitigation areas for oil and gas vessels."d <-read_csv(risk_reduction_by_areas_csv) |>select(yr, matches("wan|wab$")) |>rename_with(~str_replace(.x, "^pct_risk_reduced", "pctriskreduced"),starts_with("pct_risk_reduced")) |>pivot_longer(-yr,names_to ="var",values_to ="val") |># table(d$var)separate_wider_delim( var,delim ="_",names =c("v","type","speed","area")) |>pivot_wider(names_from =c(v, area),names_sep ="_",values_from = val) |>relocate(speed, type) |>arrange(speed, type, yr)# names(d) |> sort() |> paste(collapse="\n") |> cat()d <- d |>mutate(yr =as.character(yr)) |>bind_rows(bind_cols(yr ="Avg", d |>select(-yr) |>group_by(speed, type) |>summarise(across(where(is.numeric), mean),.groups ="drop") ) ) |>arrange(speed, type, yr)d <- d |>mutate(speed =case_match( speed, "gt01"~"All speeds","gt10"~"> 10 knots"),type =case_match( type, "all"~"All vessels","boem"~"Oil & Gas vessels") )d |>group_by(speed, type) |>gt() |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent(c( pctriskreduced_wab, pctriskreduced_wan), decimals =0) |>opt_stylize(style =6) |>cols_label(yr ="Year",risk_wab ="Original",pctriskreduced_wab ="%",risk_wan ="New",pctriskreduced_wan ="%") |>tab_spanner(label ="Risk Reduction by Area",columns =c( risk_wab, pctriskreduced_wab, risk_wan, pctriskreduced_wan))```## Reproducible Results## References {.unnumbered}